Coef-plot

library(tidyverse)
library(rstan)
library(bayesplot)

mcmc_intervals2 <- function(data,
                           transformations = list(),
                           ...,
                           prob = 0.5,
                           prob_outer = 0.9,
                           point_est = c("median", "mean", "none"),
                           rhat = numeric()) {
#  check_ignored_arguments(...)


  color_by_rhat <- rlang::has_name(data, "rhat_rating")
  no_point_est <- all(data$point_est == "none")

  x_lim <- range(c(data$ll, data$hh))
  x_range <- diff(x_lim)
  x_lim[1] <- x_lim[1] - 0.05 * x_range
  x_lim[2] <- x_lim[2] + 0.05 * x_range

  # faint vertical line at zero if zero is within x_lim
  layer_vertical_line <- if (0 > x_lim[1] && 0 < x_lim[2]) {
    vline_0(color = "gray90", size = 0.5)
  } else {
    geom_ignore()
  }

  args_outer <- list(
    mapping = aes_(x = ~ ll, xend = ~ hh, y = ~ parameter, yend = ~ parameter),
    color = get_color("mid")
  )
  args_inner <- list(
    mapping = aes_(x = ~ l, xend = ~ h, y = ~ parameter, yend = ~ parameter),
    size = 2,
    show.legend = FALSE
  )
  args_point <- list(
    mapping = aes_(x = ~ m, y = ~ parameter),
    data = data,
    size = 4,
    shape = 21
  )

  if (color_by_rhat) {
    args_inner$mapping <- args_inner$mapping %>%
      modify_aes_(color = ~ rhat_rating)
    args_point$mapping <- args_point$mapping %>%
      modify_aes_(color = ~ rhat_rating,
                  fill = ~ rhat_rating)
  } else {
    args_inner$color <- get_color("dark")
    args_point$color <- get_color("dark_highlight")
    args_point$fill <- get_color("light")
  }

  point_func <- if (no_point_est) geom_ignore else geom_point

  layer_outer <- do.call(geom_segment, args_outer)
  layer_inner <- do.call(geom_segment, args_inner)
  layer_point <- do.call(point_func, args_point)

  # Do something or add an invisible layer
  if (color_by_rhat) {
    scale_color <- scale_color_diagnostic("rhat")
    scale_fill <- scale_fill_diagnostic("rhat")
  } else {
    scale_color <- geom_ignore()
    scale_fill <- geom_ignore()
  }

  ggplot(data) +
    layer_vertical_line +
    layer_outer +
    layer_inner +
    layer_point +
    scale_color +
    scale_fill +
    scale_y_discrete(limits = unique(rev(data$parameter))) +
    xlim(x_lim) +
    bayesplot_theme_get() +
    legend_move(ifelse(color_by_rhat, "top", "none")) +
    yaxis_text(face = "bold") +
    yaxis_title(FALSE) +
    yaxis_ticks(size = 1) +
    xaxis_title(FALSE)
}

#setwd("~/seedling-stan")
files <- list.files("../rda")
#files <- files[1:15]

1 loo

n_samp <- NULL
#N <- 2
for (n in 1:length(files)) {
#for (n in 1:N) {
  file_name <- str_c("../rda/", files[n])
  file_name
  load(file_name)
  assign(files[n], fit)
  n_samp[n] <- list_dat_d$N
}

dat <- tibble(fit = files) %>%
# head(2) %>%
 mutate(N = n_samp) %>%
 mutate(elpd_list = map(fit, ~ loo(get(.)))) %>%
 mutate(looic = map_dbl(elpd_list, ~ .$estimates[3, 1]))

rm(list = ls()[ls() %>% str_detect("\\.rda$")])

DT::datatable(dat)

2 diagnostic

• skip

color_scheme_set("mix-brightblue-gray")
check_trace <- function(x) {
  np_cp <- nuts_params(fit)
  mcmc_trace(fit,
             pars = c("lp__"),
             regex_pars = "gamma",
             np = np_cp) 

  ggtitle(obs_files[x])
  }
}
p <- 2  %>%
  map(check_trace)
p %>%
  walk(print)

3 coef-plot

season-habitat-trait_data

files <- list.files("../rda")
#N <- 2
for (n in 1:length(files)) {
#for (n in 1:N) {
  file_name <- str_c("../rda/", files[n])
  file_name
  load(file_name)

  tmp_name <- trait6 %>%
    select(-sp) %>%
    names 

  tr_name <- c("Int", tmp_name)
 
  title <- str_c(dry, "-", hab, "-", trait_data)

  x_name <- names(Xd)

  gamma_name <- str_c(
    rep(tr_name, length(x_name)),
    rep("-", length(tr_name) * length(x_name)),
    rep(x_name, each = length(tr_name)))

  fig_dat <- mcmc_intervals_data(
    fit,
    regex_pars = "gamma",
    #prob = 0.5,
    prob_outer = 0.95,
    point_est = "median"
  ) %>%
    mutate(parameter2 = gamma_name) %>%
    filter(parameter2 != "Int-Int") %>%
    mutate(sig = ifelse(ll * hh < 0, "non-sig", "sig"))

  fig_dat

  p <- ggplot(fig_dat) +
    geom_vline(xintercept = 0) +
    geom_segment(data = fig_dat,
      aes(x = ll, xend = hh, y = parameter2, yend = parameter2),
      col = "darkgrey") +
    geom_point(data = fig_dat,
      aes(x = m, y = parameter2, col = sig),
      size = 4) +
    xlab("Coef") +
    ylab("") +
    ggtitle(title) +
    theme_classic()
  print(p)

}